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We show that two photons coupled to Rydberg states via electromagnetically induced trans¬ 
parency can interact via an effective Coulomb potential. This interaction gives rise to a continuum 
of two-body bound states. Within the continuum, metastable bound states are distinguished in 
analogy with quasi-bound states tunneling through a potential barrier. We find multiple branches 
of metastable bound states whose energy spectrum is governed by the Coulomb potential, thus 
obtaining a photonic analogue of the hydrogen atom. Under certain conditions, the wavefunction 
resembles that of a diatomic molecule in which the two polaritons are separated by a finite “bond 
length.” These states propagate with a negative group velocity in the medium, allowing for a simple 
preparation and detection scheme, before they slowly decay to pairs of bound Rydberg atoms. 

PACS numbers: 42.50.Nn, 32.80.Ee, 34.20.Cf, 42.50.Gy 


Photons are fundamental massless particles which are 
essentially non-interacting for optical frequencies. How¬ 
ever, a medium that couples light to its atomic con¬ 
stituents can induce interactions between photons. This 
interaction may lead to exotic, many-body states of light 
W, or can be used as a basis for realizing deterministic 
quantum gates between two photons Hd. A promising 
approach to create strongly interacting photons is to cou¬ 
ple the light to atomic Rydberg states 13 H EH HHM! , as 
realized in recent experiments 37 52] . 

Rydberg polaritons are superpositions of Rydberg 
atoms and light, which propagate almost without dissipa¬ 
tion under the conditions of electromagnetically induced 
transparency (EIT) 8. f53l [55] . EIT strongly reduces the 
group velocity and makes Rydberg polaritons dispersive. 
The large admixture of the Rydberg state can induce 
strong interactions between polaritons via the inherent 
Rydberg-Rydberg interactions. Specifically, the block¬ 
ade effect prevents the formation of two Rydberg polari¬ 
tons within the so-called “blockade radius” of each other 
[38] H3l 1481 I56H59] . When the probe photons are detuned 
from the atomic transition, they can form bound states. 
A shallow bound state of light was observed in recent ex¬ 
periments [45], while stronger interactions result in deep 
bound states of Rydberg polaritons tied together within 
the blockaded region [29] . One can imagine these bound 
states as consisting of a photon trapped by a Rydberg 
excitation in a deep square well. 

In this Letter, we predict and explore a class of 
photonic states resembling diatomic molecular states in 
which the two bound photons can be separated by a non¬ 
zero “bond length.” This is achieved by considering Ry¬ 
dberg polaritons with the quantized light red-detuned 
from the excited atomic state. In such a system, we 
show the existence of metastable bound states exhibit- 



Figure 1: (a) The probe field couples the ground state | g) 
to the excited state |e) and is red-detuned by A. A con¬ 
trol held with Rabi frequency H couples |e) to the Rydberg 
state | r) and is blue-detuned by A, thus putting the probe on 
an EIT transmission resonance. The Rydberg state is thus 
shifted downward by H 2 /A. The van der Waals interaction 
with another reference Rydberg excitation at r = 0 can bring 
| r) into an absorption resonance with the two-photon transi¬ 
tion. (b) The effective potential of two Rydberg polaritons 
as a function of their separation r. At large separations, the 
effective potential is that of the van der Waals interaction; 
the resonant condition near the blockade radius gives rise to 
a singularity, while at small separations the interaction levels 
off as the Rydberg state is highly shifted out of resonance. In 
the vicinity of the singularity, the potential behaves as that 
of the Coulomb interaction. 


ing the Coulomb spectrum, akin to the hydrogen atom. 
Such states can potentially be used as building blocks for 
more complex quantum states of light. 

To gain an intuitive understanding, consider the level 
structure of the Rydberg medium shown in Fig. 0a). 
The probe field coupling the ground state | g) to the in¬ 
termediate excited state \e) is red-detuned by A > 0, and 
the Rabi frequency of the control field coupling \e) to the 
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Rydberg state | r) is U. For H<A, the Rydberg state is 
shifted downward by U 2 /A [see Fig. [lja)]. The van der 
Waals interaction V(r) = Cq/v 6 between Rydberg states 
modifies this picture (we assume Cq > 0 or more gen¬ 
erally Cq A >0). In particular, at small separations r, 
the strong interaction shifts two Rydberg states upward 
and out of resonance, while at large separations, the in¬ 
teraction is negligible and the energy level of each atom 
asymptotes to —U 2 /A (we set h = 1). For intermediate 
separations on the order of the blockade radius 77 , de¬ 
fined by V(rb) = 2£l 2 /A (65], the system goes through 
a resonance (the factor of two arises since both atoms 
experience the U 2 /A shift). This resonance, associated 
with a pair (or “molecule”) of Rydberg atoms, endows 
the effective interaction V e ^(r) between two Rydberg po- 
laritons with a singularity separating repulsion outside 
the blockade region from attraction inside; see Fig. 0 b). 
This effective interaction between two Rydberg polari- 
tons can be roughly thought of as the difference in sus¬ 
ceptibility of a single Rydberg polariton with and with¬ 
out a Rydberg excitation at r = 0 [45] . Interestingly, the 
effective potential near the resonant edge is that of the 
Coulomb interaction that changes sign across the block¬ 
ade radius. This potential admits a continuum of states 
consisting of pairs of bound Rydberg atoms, i.e., Rydberg 
molecules. Within the continuum, we identify multi¬ 
ple branches of metastable states whose lifetime diverges 
with the strength of the interaction. When the effective 
energy of the two-polariton state lies below both Wff(°o) 
and Wff( 0 ), the bound state experiences a repulsive core 
and the wavefunction becomes double peaked near ± 77 , 
resembling a diatomic molecular state. We further show 
that the group velocity of these states is negative, con¬ 
sistent with the fact that they have a finite lifetime. 

Model. —To describe a propagating polariton in a Ry¬ 
dberg medium, we define <G(z) and (z) as creation 
operators for a photon and a Rydberg excitation, respec¬ 
tively, at position z. We define g to be the collectively 
enhanced atom-photon coupling |53 and assume that the 
decay rates 7 of the excited state (satisfying 7 < A) and 
y of the Rydberg state can be neglected. In the regime 
of slow light (g U) and with large single-photon detun¬ 
ing (A ^ O), one can adiabatically eliminate the excited 
state \e) Ha Eg. The two-state Hamiltonian of the Ry¬ 
dberg medium is then 


H = 


-I 


dz 


-icd z + g 2 /A Qg/A \ 
Vtg/A Q 2 /A ) \ S 


^ J dzdz'V(z — z')S\z)S\z , )S(z , )S(z). (1) 


in Eq. 0 can be projected onto the sector containing 
two-particles (at positions z and z') described by the 
quantum state |<F) with amplitudes EE(z,z '), ES(z,z'), 
SE(z,z'), and SS(z,z'), defined in terms of the vacuum 
|0) as SE(z,z') = (0|S'(z)5(z')|4>) and similarly for the 
other amplitudes. The problem is simplified by noting 
that, for two particles, the total energy uo and the center 
of mass momentum K are good quantum numbers. 

In the limit g —>> 0, the SS component decouples from 
the photonic amplitudes [ujSS(z, z') = (—2U 2 /A + U(z — 
z'))SS(z — z')\ giving rise to a continuum of (delta- 
function) states of Rydberg molecules. (We are mak¬ 
ing the frozen gas approximation, where we neglect the 
atomic motion as it is much slower than the polariton 
dynamics.) Upon increasing g the continuum of states is 
still present while the wavefunction remains localized to 
the blockade radius. To see this, note that the Heisen¬ 
berg equations of motion for the above amplitudes im¬ 
mediately lead mm to the Shrodinger-like equation 



G 

r 6 — [ r b{u)] 6 + i0+ 


ip(r) = F/0( r )> 


( 2 ) 


where r is the relative coordinate of the two particles, 
and -0 is the symmetrized light-Rydberg wavefunction 
ijj(r) = [ES(r ) + SE(r)\/ 2. Notice that the van der 
Waals potential is replaced by an effective potential 
V e s(r) = C 6 /(r 6 — [77 (a;)] 6 + i0 + ) modified within the 
blockaded region as in Fig. [TJb). For a nonzero cj, the 
blockade radius 77 ( 0 ;) depends on frequency via the reso¬ 
nant condition G/[?7G)] 6 = 2£l 2 / A-\-uo (the dependence 
of 77 on uj will often be implicit below). z0 + in V e s is ob¬ 
tained in the limit of vanishingly small 7 and 7', which is 
further required by causality. In the limit of small energy 
and momentum, m is the mass of a single dark-state po¬ 
lariton due to the curvature of linear susceptibility and 
is given by m = g 4 /2c 2 U 2 A [61] [62], while the energy is 
given by E = uj — v g K with v g = (U 2 / g 2 )c being the EIT 
group velocity. More generally, the parameters in Eq. 
0 can be simply derived from single polariton physics: 
For two Rydberg dark-state polaritons with momenta k\ 
and &2 and dispersion = ^(^ 1 , 2 ); the constraints 
uji T 0 J 2 — bj and k\ + = K yield an expression for the 

relative momentum p = \JmE consistent with the full, 
yet complicated, expressions for m and E 29 ; 6(1. It is 
worth pointing out that scattering theory techniques can 
be used to derive Eq. (pi) without adiabatically eliminat¬ 
ing the excited state [29], and to show its validity for a 
wide range of parameters [66] 


In the absence of interactions, H diagonalizes into dark- 
and bright-state polaritons, where, at low energies, the 
former (oc gS^ — when d z = 0) is mostly com¬ 
posed of | r) and travels at a reduced group velocity 
m- In the presence of interactions, the Hamiltonian 


Coulomb states .—The effective potential U e ff diverges 
as l/(r± 77 ,) near the blockade radius, like a Coulomb po¬ 
tential. Across the singularity, the wavefunction -0 should 
be continuous, while its derivative does not have to be. 
The full wavefunction = (EE, ES, SE, SS) has 
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various components which are related to i/j as 
2gQ/A 


EE(r) = - 


2g 2 / A + (jj — cK 


4>(r), 


(3) 


ES(r ) = 1 - 


SSlr) = ^V 


{g 2 + n 2 )/A+ui-cK/2 
ip(r) 


d r ) ip(r), 


r .—6 _ 

T 'b J 


+ a5(r± r&), 



where, for states with ip(r b ) ^ 0 , the principal value sym¬ 
bol V removes the 1 /[r±r b (uj)\ singularity in SS near the 
blockade radius. The coefficient a of the delta-function 
is determined by the discontinuity in the derivative of xjj 
at the blockade radius [60.. 

We now notice that Eqs. (2j3) admit a special set of 
solutions, which are a superposition of a normalizable 
wavefunction vanishing for \r\ > r b [ip(r b ) = 0 ] and a 
delta function singularity in the SS component, but with¬ 
out the l/[r ± r b (oo)] singularity. Such states can be in¬ 
terpreted in analogy to a leaky box where a quasi-bound 
particle tunnels through a potential barrier: for the leaky 
box, a true eigenstate is a superposition of the metastable 
bound state and a plane wave, which is a momentum 
eigenstate selected from a continuum. Similarly, for the 
above special eigenstates [with ip(r > r b ) = 0 ], the role of 
the continuum of eigenstates is played by the delta func¬ 
tions in SS , which are position eigenstates. When the 
delta function is removed, the other components of the 
wavefunction form a metastable bound state. Further¬ 
more, in the limit of an infinitely strong interaction, i.e. 
g oo, our special eigenstates lose their delta function 
component m- This is again analogous to the leaky box, 
where, in the limit of an infinitely tall barrier (i.e. the no¬ 
leak limit), one obtains exact eigenstates confined to the 
box and decoupled from the plane-wave component sit¬ 
ting outside the box. Henceforth, we call the metastable 
bound states above (without the delta function) Coulomb 
states, and study their spectrum and other properties in 
detail. 

Figure [2j[a) shows the energy spectrum of the exact 
eigenstates (i.e. with the delta function) underlying the 
Coulomb states. The exact solutions are depicted as solid 
lines, while the dashed lines show the energy spectrum 
derived from the WKB quantization condition [applied 
to the case ^{±.r b ) = 0 ] 



dr = n 7 r, 


n= 1 , 2 , ••• 


(4) 


with p(r) = yjm(E— V e ^ (r)) defined by Eq. ([ 2 ]), and 
r 0 < r b is the classical turning point near the origin 
EH- Figure [2j[a) demonstrates that the WKB quanti¬ 
zation agrees with the full solution for values of K near 
2g 2 /cA. 

When K is close to 2 g 2 /cA, we can analytically com- 


Figure 2: (a) Dispersion curves for the exact eigenstates un¬ 
derlying the metastable Coulomb bound states (only the first 
four branches n = 1-4 are shown) with g 2 rb/cA — 40 and 
£t/g = 0.05. The solid lines give the exact solution, while the 
dashed lines represent the WKB results. For K —>• 2g 2 /cA, 
the WKB results are almost exact. The dispersion curves 
converge to one point with a negative slope, or group veloc¬ 
ity. (b) Decomposition of the metastable Coulomb bound 
states (T$b defined as ^ aJn) x n with the delta-function con¬ 
tribution removed) into the continuum of exact eigenstates. 
Here we took g 2 rb/cA = 40 and Q/g = 0.01 with a fixed cen¬ 
ter of mass momentum cKA/2g 2 = 0.95. The width of these 
distributions is much less than the energy spacing, a strong 
indication of metastability. The inset shows the wavefunction 
components for the n — 1 Coulomb state with the param¬ 
eters in (a) and cKA/2g 2 = 0.98. (The EE component is 
exaggerated by a factor of 1.5 for better visibility.) 


pute the integral in Eq. © to find 

l+uA/20 2 [g 2 r b (w)/cA] 2 

1 - cKA/2g 2 n 2 ’ 1 j 

where A = [r(2/3)/r(l/6)y / 7r] 2 ~ 0.014. If r b were in¬ 
dependent of cj, Eq. © would imply that uj is quan¬ 
tized as 1/n 2 (plus a constant), reminiscent of the en¬ 
ergy spectrum of the Coulomb potential. However, due 
to the l o dependence of r&, the quantization changes to 
c o n ~ 1/n 3 / 2 [ 68 ]. The fact that the blockade radius, 
and, thus, the interaction strength, is sensitive to fre¬ 
quency is a typical feature of nonlinear optical systems 
m ■ We also stress that Ag 2 r b /cA is identical to the fig¬ 
ure of merit in the far-detuned regime OD b ^/ A, where 
OD b is the optical depth per blockade radius. The fig¬ 
ure of merit quantifies the strength of the interaction as 
two polaritons imprint a phase ~ OD b 7 /A on each other 
[45]. 

With the dispersion in hand, we now explore the sta¬ 
bility of the Coulomb states. The solutions given by 
Eq. © are a complete set of eigenstates for the two- 
particle Hilbert space. To normalize these states, we 
take K to be fixed and use the energy normalization 
(^cj',k\^cj,k) = S(cj — c */). We can then verify the 
metastability of the Coulomb states (^[r&(u;)] = 0 ) with 
the delta-function removed by looking at their decompo¬ 
sition into the normalized eigenstates. Figure 2(b) shows 
this decomposition for several n, where we see that the 
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Coulomb states are sharply peaked at the expected fre¬ 
quencies. The width of these distributions can be much 
narrower than the spacing between states, a strong sig¬ 
nature of spectral distinguishability [63] . Furthermore, 
the Coulomb states converge to the exact eigenstates for 
a very strong interaction strength, which is analogous to 
the leaky box in the limit of an infinitely deep potential 

EOJ- 

A unique feature of the dispersion curves in Fig. Ha) 
is that their slope, and thus the group velocity, is 
negative. While true eigenstates cannot have a nega¬ 
tive group velocity in the absence of left-moving modes 
(Supplemental Material [60]), Coulomb states are not 
exact eigenstates and eventually decay into Rydberg 
molecules. Equation © gives the group velocity as 
v = —A (# 2 n>/cA) Vg/n 2 , where v g is the EIT group ve¬ 
locity. Therefore, the velocity is also quantized as 1/n 2 
for different branches of bound states (and fixed values 
of uj). This quantization and the negative sign make the 
group velocity an ideal signature for detecting different 
Coulomb states. We also remark that a small 7 (< A) 
modifies the energy only slightly, proportional to 7/ A, 
which thus becomes negligible for large detuning. 

We now show how to prepare these states and measure 
their dispersion. We assume that we have access to an 
additional hyperfine ground state | q), which, as shown 
in Fig. [3p)> is connected to both | g) and the Rydberg 
state | r) through two-photon transitions via an excited 
state \e'). With these additional states, we can effectively 
turn on and off the polariton interactions by applying a 
fast 7r-pulse on the two-photon transition between | q) and 
|r>. 

The preparation procedure is as follows. First, we store 
two identical photons (equivalently a weak coherent state 
followed by postselection) in the atomic state | q) using 
standard protocols [53, 64]. To have a significant over¬ 
lap with the Coulomb states once we map to |r), the 
state has to have the correct center of mass momentum 
K. To achieve this, we introduce a linear energy gradi¬ 
ent E' along the atomic cloud for a time r, which could 
be achieved with a magnetic field gradient, another op¬ 
tical beam, or microwave field. This will impart a phase 
e -iE tR on s tored two-photon state. By choosing the 
appropriate r and then mapping | q) to |r), we can selec¬ 
tively excite different Coulomb states provided they have 
a large enough spatial overlap with the initial product 
state input. As the bound states travel with negative 
group velocity, the Coulomb state component will sep¬ 
arate from the rest of the wavefunction. To detect the 
state, one can then map the Rydberg state back to | q) 
and either measure the population of the state | q) di¬ 
rectly or retrieve the state into light. In Fig. [3][b-d), for 
realistic parameters (including 7 and 7') [411, we verify 
this approach by preparing a variational state that has 
a large overlap with the SS component of the Coulomb 
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Figure 3: (a) Level structure used to prepare initial SS dis¬ 
tribution. (b-d) Time evolution of a wavepacket with all 
components initially zero except SS , which is chosen to be 
a Gaussian wavepacket of variational n — 1 Coulomb state 
solutions (with the delta function removed) centered at uj m 0 
and having width Q 2 /2A ^0i. Specifically, \EE\, initially 
zero, is shown after the initial transient evolution subsides 
at (b) t = tf/ 4, and at (c) t = tf/2 and (d) t — tf , where 
tf = 10 A/O 2 . The wavepacket within the blockade radius has 
the expected shape of the Coulomb state, propagates back¬ 
ward, and decays, while the wavepacket outside the blockade 
radius propagates forward with v g and disperses. We took 
a medium of length L — 1 677 , g 2 r b /cA = 5, g/2n = 17 
GHz, 0/2 tt = 1.5 MHz, r b = 25 /im, A/2 t r = 30 MHz, 
7 / 27 r — 3 MHz, and 7 / /27 t — 100 kHz [451 60 . 


state [shown in Fig. [3](b)] with other components equal 
to zero and solve numerically for the time evolution 60 . 
In this case, the effective energy E of the bound state lies 
above 4/^(0) and the wavefunction is peaked at r = 0. 
We have also verified that, when E < (0), the back¬ 

ward propagating state becomes double peaked and that, 
for smaller decay rates and larger g 2 77 /cA, the negative 
group velocity observed in the numerics agrees with the 
WKB prediction from Eq. © for the n = 1, 2 and 3 
Coulomb states within a few percent in each case m- 
Before concluding, we point out several conditions 
necessary for the experimental observation of Coulomb 
states. First, the delta functions appearing in the un¬ 
derlying exact eigenstates should be consistent with the 
treatment of atoms as a continuous medium. This re¬ 
quires p(ttw 2 )ws 1, where p is the atomic density, 
w is the beam waist, and w$ ~ \dr b (uj)/duj\/r is the 
effective width of the delta function due to the uncer¬ 
tainty in uj coming from the finite duration of an ex¬ 
periment r, which is in turn limited by the lifetime of 
the Coulomb states A/2f) 2 from Fig. [ 2 ) 3 ]. Recent 
experiments achieved a density ofp = 2xl0 12 cm -3 
with a beam waist of w = 4.5 /im, corresponding to 
g/27r ~ 4 GHz [45]. We will use below the parameters of 
Fig. i where g/2'K = 17 GHz, so we will take w = 4.5 pm 
and p = 3.6 x 10 13 cm -3 . Taking r = £/, we then find 
ws ~ 0.2 pm and p(ttw 2 )ws ~ 200. Second, two atoms 
initially 77 away from each other should not change their 
distance by more than w$ if they are displaced trans¬ 
versely by w. This leads to the condition w < 
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which is nearly satisfied, and can be more strictly satis¬ 
fied by changing the center frequency of the wavepacket 
to increase 77 ,( 0 ;). Finally, the force on a pair of Ry¬ 
dberg atoms r\) apart and their thermal velocity must 
both induce motion less than ws over time r, leading 
to the conditions 6Cqt 2 / (mrl), rksT/m < w§. Using 
the mass m of 87 Rb and temperature T — 35 fiK [41], 
the first condition is satisfied (0.07 fim < 0.2 fim). The 
second condition can be satisfied by using a sufficiently 
high control field intensity. 

Outlook .—While our proposal opens the avenue for the 
creation of Coulomb-like two-photon states, we expect 
that a wide class of both useful and exotic two-photon 
and multi-photon states can be created via refined engi¬ 
neering of photon-photon interactions, e.g. by using mi¬ 
crowave fields [43]. The detailed understanding of the 
two-photon Rydberg-EIT physics provided by this work 
also opens up an avenue towards understanding the full— 
and much richer—many-body problem involving an ar¬ 
bitrary number of photons in any dimension. 
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Supplemental Material for “Coulomb bound states of strongly 

interacting photons” 


I. ADIABATIC ELIMINATION 

In this section, we derive the Shrodinger equation (2) 
of the main text. 

Upon adiabatic elimination of the intermediate state 
|e), the two-particle wave function is described by four 
components EE, ES, SE , and SS , each of which is a 
function of time and two spatial coordinates. We define 
ES± = (ES d= SE)/ 2 as the symmetric and antisymmet¬ 
ric combinations of the photon-Rydberg components. In 
the frequency domain, the Heisenberg equations can be 
cast as (see Refs. [Sl ] -[S4 ] for more details) 


which, for CqA > 0, reduces to the effective potential 
in Eq. (2) of the main text. Defining the normalized 
units O = cdA/2D 2 and K = cKA/2g 2 , the values of the 
energy and the mass take the form 


„ 2D 2 _ n2 

E — -^-(1 + (j) 


Q 2 

1 — K -\ - — (\ + 20) — — 

9 1 


n 2 /g 2 


(S10) 

1 


K + on 2 /g 2 


1 


m = 


2D 2 Ac 2 (1+cd) 2 


Q 2 

1 — K -\ -^-(1 + 2cu) 

9 


1 + UJ 

(Sll) 


to EE = —icd R EE -E-EE- ^ ES +, 

luES+ = ~d R ES+ - E^ ES+ 

- icd r ES _ - ^ (EE + SS), 


uiES- 


- l /d R ES _ 


g 2 + n 2 


ES- - icdrES. 


+> 


2Q 2 2qQ 

uSS = ——SS + V(r)SS --j^ES+, 


(51) 

(52) 

(53) 

(54) 


where uj is the frequency, r = z — z' denotes the relative 
coordinate, and R = (z + z')/2 is the center of mass co¬ 
ordinate. For an infinitely long medium, one can work 
in Fourier space (relative to R) with the total momen¬ 
tum K. Defining ^p(r) = ES+(r ), Eqs. flSl|S3|S4D yield, 
respectively, 


EE = - 


ES- = 


2 gtt 


1 


A 

A 


uj — cK 




—ic 


. CJ _ Kc 

A \ dJ 2 


drip, 


„„ 2 gfi 

SS = — p—v 

A 




2 ~¥+oo-V(r) 


(55) 

(56) 

+a5[r - r h (uj)\, (S7) 


with V denoting the principal part near the singularity 
at the blockade radius. The coefficient a is determined 
by matching boundary conditions across the singularity. 
Inserting these expressions into Eq. (S2), we obtain a 
second-order differential equation for ^ as 


- -d 2 ^ + V eff (r)V> = Exp, 

m 


(S8) 


which is valid everywhere away from the blockade radius, 
i.e. for \r\ > r^uj) and \r\ < r^uj). The effective poten¬ 
tial is given by 


V etf (r) 


V(r) 


l-V(r)/( 2 ^+u)' 


(S9) 


These expressions also fully agree with the diagrammatic 
approach of Ref. [S4] in the above limits. 


The boundary conditions at the origin [^'(O) = 0] and 
at infinity [ip(r —>> oc) = 0] are necessary to solve for 
%b. Furthermore, the wavefunction should be continuous 
across the singularity at r = r&. On the other hand, the 
discontinuity in its first derivative at r — determines 
the coefficient a , via Eqs. flS2|S3|S7D , as 


Ac/gVt 

(g 2 + Q?)/Ac + oj/c — K/2 



(S12) 


II. BOUNDARY CONDITION AT THE 
SINGULARITY 


In this section, we show how to ex plici tly calculate 
d r 'ip \ r -, and hence determine a via Eq. (S12), taking into 

r b I 

account the boundary conditions at r = U and r = oo. 
This is necessary for constructing the eigenbasis used in 
Fig. 2(b) of the main text. 

Because of the singularity in Wff, we find that a has 
both real and imaginary parts which can never simul¬ 
taneously vanish. This implies that all eigenstates have 
a delta-function contribution. To show this, consider a 
small neighborhood near the singularity where Eq. (S8) 
takes the form 


1 

U 2 


%'!> = 


l 


X — 1 




(S13) 


with U ~ 9 -/^E\/t^k/S- and x = r /n(u). This 
equation is valid for \x — 1| <C 1 and has analytic solu¬ 
tions on both sides of the singularity in terms of first 




























order Bessel functions 


(S14) 

^(x)^- 7 ^^Y 1 (2UVl^c), (S15) 

'iP+(x)*^^I 1 (2UV^i), (S16) 

i>t ( x ) ~ 2 ^~ 1 K\(2U\fx - 1), (S17) 


where =p refers to x ^ 1, and the ~ signs are meant to 
indicate that these equalities hold only for \x — 1| <C 1. 
The solutions ip for all x are obtained by using Eq. (S8) 
to extend the above ip^ 2 t° x - One then imposes the 
boundary conditions that 'ip is symmetric about x = 0, 
i.e., dip/dx(x = 0) = 0, and fp(x) —» 0 as x —>• oo. The 
solution for x < 1 can be written as 


ip{x) = c 1 ip 1 (x) + c 2 tp 2 (x), (S18) 


where c\ and C2 are determined by the boundary condi¬ 
tion at x — 0. For x > 1, one can show that the choice 
of ^ i n terms of the Bessel function of the second kind 
K\ near the singularity guarantees that, away from the 
singularity where Eqs. (S15)-(S17) are no longer valid, 
ip2 decays exponentially for large values of x, while 
grows exponentially. Imposing the boundary condition 
that 'ip vanishes at infinity implies ip(x) oc ip£(x), while 
the continuity of ip(x) at x = 1 gives the result 


^(x) 


CiVq (x) + C 2^2 ( X )i X < 1? 

02^2 ( X )i X > 1 . 


(S19) 


The contribution of ipf to d x ip\\- is c i- However, cal¬ 
culating the contribution of 'ipf to d x ip\\+ requires more 
care because, near the singularity, d x ip± ~ log(l — x). To 
help resolve this, we examine the solutions in the presence 
of an infinitesimal decay rate 7' from the Rydberg state. 
In this case, near the singularity, the effective Schrodinger 
equation takes the form 

h 3 ^ = (S20) 

where e = 7 / /( Co ’ + 2U 2 /A). Using the relation 

lim-= V\ - ] — iirStx — 1), (S21) 

e-^o x — 1 T ie \x — l) 

and integrating Eq. ( |S20| ) across the singularity suggests 
a contribution to d x ip\\_ equal to —i 7 rU 2 r ip(l) = —11:02. 
Combining this result with the c\ contribution gives 


a = — 


A c/gVtr b 


( 1 g 2 + Q 2 )/Ac + uj/c — K/2 


(ci—i7TC2). (S22) 


With this final result, we can construct the complete set 
of states used in making Fig. 2(b) of the main text. First, 


we solve numerically for ip inside the blockade radius. 
We then find c\ and C2 by matching these numerical so¬ 
lutions, in a region near the singularity, to the analytic 
solutions in Eqs. (S14)-(S15). To check the arguments 
presented above, we have also verified numerically that 
when the imaginary component to a is neglected the re¬ 
sulting solutions do not form an orthonormal basis, while, 
when the imaginary term is included, the solutions are 
consistent with an orthonormal basis. 

Since can be assumed to be real numbers, which 
cannot both vanish, it is clear that all eigenstates have 
the delta-function contribution. This is consistent with 
the fact that, for g = 0, the continuum is composed of 
Rydberg molecule states. For finite g , these states be¬ 
come dressed with the photons, but do not lose their 
character as atomic bound states. The eigenstates linked 
to the Coulomb states have the special property that 
c 2 = U 2 ^P( 1) = 0. 


From this solution, we can also determine the behav¬ 
ior of the eigenstates as the interaction strength U 00. 
At the singularity, = C2/U 2 , which implies that, as 
U increases, ip( 1) —» 0 and all solutions satisfy the same 
boundary condition at the singularity. Additionally, in 
this limit a 0, which implies that the states confined 
inside the Rydberg blockade completely decouple from 
the continuum of Rydberg molecule states (i.e., the delta 
functions). This is again analogous to the leaky box dis¬ 
cussed in the main text, where, as the box becomes in¬ 
finitely deep, the eigenstates inside the well become de¬ 
coupled from the continuum of momentum states that 
live outside the box. 


III. NUMERICAL METHODS 

In this section, we describe the numerical methods used 
to obtain Fig. 3 of the main text. 

We include the decay rate 7 of the intermediate 
state by adding an imaginary component to A. De¬ 
cay of the Rydberg state requires adding the term 
—i{p( , /2) f dz S 1f (z)S(z) to the Hamiltonian. 

Within the two-excitation subspace, H can be split 
into a kinetic term T that describes the propagation of 
photons and a part W that is diagonalizable in real space, 
that includes the Rydberg-Rydberg interaction, decay, 
and coupling to the quantum and classical light fields. 
We can then find the time-evolution of the wavefunction 
by a Trotter decomposition, whereby we split the prop¬ 
agator into two parts which are separably diagonalizable 
in momentum (T) and real-space (W) 

e~ iHT « e -iTr/2 e -iWT e -iTr/2 + Q([W,T } T 2 ). (S23) 

In our case, £ has a linear dispersion, which implies that 
propagation with T corresponds to a uniform shift in real- 
space of the ^-components of the wavefunction, while 
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Figure SI: Time evolution of the n = 1 Coulomb state as in 
Fig. 3 of the main text, except here we take a much larger 
value of g 2 r b /cA and much smaller decay rates to verify that 
our analytical theory accurately describes the Coulomb states. 
The initial condition for EE, ES , and SE is chosen to be zero, 
while SS is chosen to be given by Eq. (S24) with a = Q 2 /2A 
and n — 1. The \EE\ component is shown (a) shortly after 
t — 0 at tVg/L — 10 -4 and at later times (b) tv g /L = 3 • 10 -3 
and (c) tvg/L — 6 • 10 -3 . Here L is the length of the medium, 
and we took g 2 rb/cA — 40, Q/g — 0.05, L/rb — 14, Q/A — 
0.25, 7 /A = 0.05, and 7 '/A = 0.01. For these parameters, in 
contrast to those of Fig. 3 in the main text, the EE component 
of the Coulomb state is localized near the blockade radius. 


e -iWr can k e f oun d exactly for each point in space. Using 
these solutions, we can construct the long-time-evolution 
by stepwise application of Eq. (S23) for small r. 


IV. CONDITION FOR REPULSIVE CORE 


In this section, we establish the parameter regime when 
the effective energy E < Wff(0). In this case, the two- 
photon state feels a repulsive core and becomes peaked 
near the points ±77. From Eq. ( |S10| ) and Eq. (SI 1) we 
can see that, for Q 2 / g 2 <C 1, we can rewrite 


E - U eff (0) 


2Q 2 (1 


A 1 -K 


■[(1 -Kf-tf/g 2 ], 


which becomes negative when 1 — K < Q/g. Note that, 
in order for the adiabatic elimination to be valid near 
K = 1, we also require 1 — K > U 3 /A 3 [S4]. This sets 
the constraints Q/g > 1 — K >>U 3 /A 3 , in order to have 
a repulsive core, which can be easily satisfied. 

In Fig.[Sl]we show the resulting backward propagating 
state under the same preparation procedure as described 
in the main text and Fig. 3, except in this case we took 
smaller decay rates, a larger value of g 2 r b /cA and a much 
larger control field intensity. In particular, the Coulomb 
state we prepared had 1 — K ~ 0.02 and Q/g = 0.05. The 
double peaked structure is clearly visible, consistent with 
the presence of the repulsive core for this bound state. 


V. GROUP VELOCITY OF COULOMB STATES 


For the experimental parameters in Fig. 3 of the main 
text, the group velocity v g /c ~ Q 2 /g 2 ~ 10 —8 , which 
implies that there is a large separation of time scales 
between the light propagation and the atomic dynam¬ 
ics. Just increasing the time step cannot overcome this 
because the error term in the Trotter decomposition be¬ 
comes very large. This problem can, however, be over¬ 
come by bringing the time scales closer together through 
the scaling transformation z —>> (z, g g/y/C, and 
rb C r b- One can see from Eqs. ( |S10||S11| ) that the dy¬ 
namics are invariant under this transformation provided 
v gC/ c <C (1 — K) 2 . To obtain Fig. 3 of the main text, we 
use (= 1.2 x 10 7 , which satisfies this condition. 


The observation of Coulomb states requires large in¬ 
teraction strengths (equivalently atomic densities). This 
in turn requires a fine numerical mesh for the wavefunc- 
tion, which makes the simulations very time consuming. 
However, since the Coulomb states are confined to a re¬ 
gion on the order of the blockade radius we can force the 
wavefunction to be zero outside a region on the order of 
a few blockade radii. This significantly reduces the re¬ 
quired memory and simulation time. We have verified 
that our numerical results are insensitive to this cutoff. 


Finally, H is an effective Hamiltonian obtained after 
adiabatically eliminating the intermediate state \e). We 
have also verified numerically that our results hold when 
state \e) is explicitly included in the simulations. 


In this section, we compare the group velocity pre¬ 
dicted by the WKB treatment in the main text with 
numerical simulations for n = 1, 2, and 3. To con¬ 
struct the Coulomb wavepacket, we first choose a nar¬ 
row range of frequencies around u) = 0 and, for every 
cj, we find K n (uj) from Eq. (4) in the main text, which 
gives us an expression for p n (r, uj). The initial state has 
EE = ES = SE = 0 with SS given by the variational 
wavefunction 


SS n (r,R) =M [ doje iK "^ R -“ 2 /° 2 - ) 

J 1 - [r b (uj)/r} Q 

x ( l) n cos / Pn(r,u) 0[r-r 6 (o;)], 

Jo 


(S24) 


which vanishes at r b (uj) for every u). Here AT is a normal¬ 
ization constant, a is the width of the wavepacket, and 
0 is the Heaviside step function. 

The results are shown in Fig. 3 of the main text and 
in Figs. [Sl]|S2l Fig. 3 of the main text and Fig. [Si] show 
snapshots of the backward-propagating n — 1 wavefunc¬ 
tion. While Fig. 3 of the main text uses experimentally 
realistic parameters, Fig. SI assumes larger g 2 r b /cA and 
smaller decay rates to verify that our analytical theory 
describes the Coulomb states accurately. Indeed, in Fig. 
|S2| which uses the parameters of Fig. |S1[ we see excel¬ 
lent agreement between the numerical simulations and 
the predicted group velocity. For each n, the extracted 
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Figure S2: Time evolution of the average center of mass posi¬ 
tion of the Coulomb wave packets with parameters as in Fig. 


SI The initial state is given by Eq. (S24). The horizontal 
axis is propagation time in units of L/v g , while the vertical 
axis is defined as the average over r < rb of the position of 
the peak values of the EE-component of the wavefunction, 
where, at each r, the peak value is defined with respect to 
R. The n — 2 and 3 curves are shifted vertically for visibil¬ 
ity. The dashed lines are the prediction for the group velocity 
from Eq. (5) in the main text and are also plotted with a shift 
relative to solid curves for visibility. The extracted slope from 
the linear region of the simulations agrees with the predicted 
group velocity to within a few percent for each n. 


slope in the linear region agrees with the predicted value 
to within a few percent. For n = 1, 2 and 3, the group 
velocity of the Coulomb states for these parameters is ap¬ 
proximately —50 • v g , —20 • v g and —10-v g , respectively. 

VI. NON-NEGATIVITY OF GROUP VELOCITY 

In this section, we show that the group velocity in a 
system with, possibly interacting, right-going modes can¬ 
not be negative for normalizable eigenstates. Let us write 
the Schrodinger equation as 


(S25) 


where we have made explicit the dependence on some 
parameter A, which will be identified later as the total 
momentum of two particles. One can then see that 

9rlvr = r\9rHr\^ K ) + UrOr^ r\^ r) 

m (^r\9rH k |T#), (S26) 

where the second term in the first line vanishes for a 
normalized state (^r\^r) = 1. This result is known as 
the Hellmann-Feynman theorem. 

For our system, we can cast the Hamiltonian 
in the two-particle sector in the basis defined by 
(EE EE + ES- SS) T . (The generalization to the case 
where \e) is not adiabatically eliminated is straightfor¬ 
ward.) Identifying K with the total momentum, from 
Eqs. ( |Sl]|S4| ), we have 


8rH k = 5{r - r') 


fc 0 0 0\ 

0 c/2 0 0 
0 0 0 0 ’ 
\0 0 0 0 / 


(S27) 


which is a nonnegative matrix acting on the Hilbert space 
associated with the relative coordinate r (the center- 
of-mass plane wave, although not normalizable, does 
not enter the argument). Therefore, the group velocity 
v g = dujR/dK given by Eq. (S26) cannot be negative for 
any eigenstate whose relative-coordinate wavefunction is 
normalizable. 
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